SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN 
METHOD FOR FRACTIONAL DIFFUSION AND WAVE EQUATIONS* 

KASSEM MUSTAPHAt AND WILLIAM MCLEAN* 

Abstract. We consider an initial-boundary value problem for dtu — 9(~°V 2 ii = /(t), that 
is, for a fractional diffusion (—1 < a < 0) or wave (0 < a < 1) equation. A numerical solution 
is found by applying a piecewise-linear, discontinuous Galerkin method in time combined with a 
piecewise-linear, conforming finite clement method in space. The time mesh is graded appropriately 
near t = 0, but the spatial mesh is quasiuniform. Previously, we proved that the error, measured 
in the spatial L2-norm, is of order k 2+a - + h 2 £(k), uniformly in t, where k is the maximum time 
step, h is the maximum diameter of the spatial finite elements, a_ = min(«, 0) < and £ (k) = 
max(l, | logfc|). Here, we generalize a known result for the classical heat equation (i.e., the case a = 
0) by showing that at each time level t n the solution is superconvergent with respect to fc: the 
error is of order (fc 3 "*" 2 "- + h?)i(k). Moreover, a simple postprocessing step employing Lagrange 
interpolation yields a superconvergent approximation for any t. Numerical experiments indicate that 
our theoretical error bound is pessimistic if a < 0. Ignoring logarithmic factors, we observe that the 
error in the DG solution at t = t n , and after postprocessing at all t, is of order fc 3 + a - + h 2 . 
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1. Introduction. In previous work [22j [30j [3lJ [32] , we have studied discontin- 
uous Galerkin (DG) methods for the time discretization of the abstract intial value 
problem 

u' + B a Au = f(t) forO<t<T, with u(0) = u , (1.1) 

where u' — du/dt and B a = 9 t _Q ; more precisely, letting uj^{t) = i^ -1 /T(fi) for [i > 0, 
the function B a v is either a (Riemann-Liouville) fractional order derivative in time, 

d /"* 

B a v(t) = — u 1+a (t - s)v(s) ds if-l<a<0, (1.2) 
ot Jo 

or a fractional order integral in time, 

B a v{t) = I u) a (t - s)v(s) ds if < a < 1. 
Jo 

In Section [5] we set out technical assumptions on the operator A, but for the present 
discussion we simply take Au = — V 2 u on a spatial domain fi, and impose homoge- 
neous Dirichlet boundary conditions on u. 

Problems of the form arise in a variety of physical, biological and chemical 
applications [12l [HI [26l [27l [34l [38l [39l [40]. The case -1 < a < describes slow 
or anomalous sub-diffusion and occurs, for example, in models of fractured or porous 
media, where the particle flux depends on the entire history of the density gradient Vit. 
The case < a < 1 describes wave propogation in viscoelastic materials [THl El [3S] . 
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In the limit as a 0, the evolution equation in (|l.ip becomes u' + Au = /, which 
is just the classical heat equation, and Eriksson et al. [5] studied the convergence of 
the DG solution U(t) w u(t) in this case. For a maximum time step k, and using 
discontinuous piecewise polynomials of degree at most q — 1 in t, with no spatial 
discretization, they proved an optimal convergence rate 

\\U(t) - u(t)\\ < Cfc«(||«o||, + ll« (9) (0)|| + ||/ (g - 13 (0)|| + fj \\f (q) { S )\\ ds 

for < t < T, where \\v\\ is the norm in £2^) and \\v \\ q = \\A^ 2 v\\ for v G D{Ail' 2 ). 
In addition, they proved that the DG solution is superconvergent at the nth time 
level t n , satisfying an error bound 

\\U(t-) u(t n )\\ < Ck 2 ^ 1 (jKlk-r + h (9) (0)|| g -i + fj \\f (q) { S )\\ q -i ds^j , 

where U(t~) = lim t ^ 4 - U(t) denotes the limit from the left. Ericksson et al. were also 
able to prove that a convergence rate faster then 0(k q ) holds under less restrictive 
spatial regularity requirements on the solution u. Our aim is to establish supercon- 
vergence results for the fractional-order problem restricting our attention to the 

piecewise- linear DG method (q — 2). We believe our scheme is the first to achieve 
better than second-order accuracy in time. As well as nodal superconvergence of the 
DG solution we show that a postprocessed solution is superconvergent uniformly in t. 

Many authors have studied numerical methods for In the case < a < 1, 

Sanz-Serna [36] proposed a convolution quadrature scheme, and subsequently Cuesta, 
Lubich and Palencia [H \5\ [3] developed this approach to obtain an 0{k 2 ) method 
as well as a fast implementation [37] ■ McLean and Thomee [23J combined finite 
differences and quadrature in time, with finite elements in space. 

In the case — 1 < a < 0, Langlands and Henry |13j introduced an implicit Euler 
scheme involving the Griinwald-Letnikov fractional derivative and spatial finite differ- 
ences with step size h, and observed 0(k 1 ^ 2 + h 2 ) convergence in the case a — —1/2. 
Yuste and Acedo [43] treated an explicit Euler scheme and showed 0(k + h 2 ) con- 
vergence. Zhuang, Liu, Anh, Turner et al. [U [UJ |3U US] developed another class 
of 0(k + h 2 ) finite difference methods, and Yuste [32] presented an 0(k 2 + h 2 ) 
method. Cui [5] and Chen et al. [I j studied 0(k + h 4 ) schemes, and Cui [7J [5] anal- 
ysed an (9(fc mm ( 1 - Q ! 2 + Q ) 4- h 4 ) ADI scheme on a rectangular spatial domain; see also 
Wang and Wang [41] and Zhang and Sun [55] . For another type of finite difference 
scheme [281 129] . the error is 0(k 2+a + h 2 ), and recently Jin et al. [TT] proved optimal 
error bounds for two semidiscrete finite element methods. Some of these works employ 
an alternative formulation of (|1.1[) using the Caputo fractional derivative. 

In practice, the higher order derivatives of u are typically singular [THIHI] as t — > 0, 
so formally high order methods Q] [2J H 13 IE CDS 133 ED EH Hfl S3 can fail to achieve 
fast convergence. We have analysed several methods that allow for the singular be- 
haviour of u by employing non- uniform time steps [21] [25] [28] [29] [32]. Another 
approach, that yields a parallel in time algorithm with spectral accuracy even for 
problems with low regularity, is to approximate u via the Laplace inversion for- 
mula [mmm]. 

To minimise the need for handling separately the cases a < and a > 0, it is 
convenient to write a + = max(a, 0) > and a_ = min(a, 0) < for the positive and 
negative parts of a, respectively. In our theory, we assume that there exist positive 
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constants M and a such that 



| Aito 1 1 + ||Au(t)|| < M and ||Ait'(t)|| +t\\Au"(t)\\ < Mf 



XT-l 



(1.3) 



as well as 



t\\A 2 u'(t)\\+t 2 \\A 2 u"(t)\\ < Mt a 



-l 



(1.4) 



for < t < T. For instance [H [H] , if / = and u € D{A 2 ) 7 then (jTSJ and (TOl) 
hold with M = C||A 2 it || and cr = 1 + a_. 

Section[5]sets out our notation and assumptions, and recalls some tools and results 
from earlier work [31] . In Section [3J we introduce the homogeneous dual problem, 



for a given terminal value zt, and represent the nodal error U(t~) — u(t n ) in terms of 
z(t) and its DG approximation Z{t). We allow a class of non-uniform meshes, specified 
in Section 01 where we prove in Theorem 14.31 that the nodal error is 0(k 3+2a ~ ). Our 
method of analysis allows us to handle the two cases — 1 < a < and < a < 1 
together, but the former presents additional technical difficulties in some places. In an 
earlier paper [30, Theorem 4.1], we estimated the nodal error for the case < a < 1 
in a different way that yields a bound of order k 2+a . (Although we claimed 0(k 3 ) 
convergence, the first line of [301 Corollary 4.2] contains an error.) 

In Section [5] we construct, via a simple interpolation scheme, a postprocessed 
solution [/" whose error is O(k 3+2oi ~ ) for all t, not just at the nodal values. Section[6] 
introduces a fully discrete scheme by applying a continuous piecewise-linear, finite 
element method for the spatial discretization. Thus, the fully discrete solution is 
continuous in space but discontinuous in time. We show that the error bound is as for 
the semidiscrete method but with an extra term of order h 2 . Finally, we present some 
numerical examples in Section[7J which indicate that our error bounds are pessimistic, 
at least in some cases. We observe that the nodal error from the time discretization 
is 0(k 3+a ~), which is better than our theoretical estimate by a factor k a ~ . The same 
is true for the postprocessed solution, uniformly in t. 

2. Preliminaries. 

2.1. Assumptions on the spatial operator. We assume as in earlier work [£l 
[31] that the self-adjoint linear operator A has a complete eigensystem in a real Hilbert 
space H, say A(f>j = Ajtfij for j = 1, 2, 3, . . . , and that A is strictly positive-definite 
with the eigenvalues ordered so that < Ai < A2 < ■ ■ • . (Strict positive definiteness 
is not essential, but allowing Ai = would result in some technical complications that 
we prefer to avoid.) We denote the inner product of u and v in H by (u,v) and the 
corresponding norm by ||it|| = y/ (u, it). Associated with the linear operator A is a 
bilinear form, denoted by the same symbol: 



z' + B* a Az = for < t < T, with z{T) = z T , 



(1.5) 



00 




m— 1 



These assumptions hold, in particular, if A = —V 2 subject to homogenous Dirichlet 
boundary conditions on a bounded domain CI, because A has a compact inverse on H = 
Li2(Cl) and A(u, v) = J Q Vit • Vv dx. 
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2.2. The discontinuous Galerkin time discrectization. Fixing a time in- 
terval [0, T], we introduce a mesh for the time discretization, 

Q = t <*! < ••• <t N = T, (2.1) 

with k n = t n — and /„ = (i n _ i,t n ) for 1 < n < N, and a maximum time 
step fc = maxi< n <jv k n . Let P r denote the space of polynomials of degree at most r 
with coefficients in D(A X / 2 ), and let J„ — [j" =1 Ij — [0,t n ] \ {t , t\, . . . , t n }, with 
J = Jn- Our trial space W consists of the piecewise- linear functions U : J — > D{A 1 / 2 ) 
with U\j n 6 Pi for 1 < n < N. We treat U as undefined at each time level t n , and 
write 

in = U(t-), Ul = U{t+), [U] n = Ul-Ul. (2.2) 

For r G {0, 1,2,.. .} we let C r (J, H) denote the space of functions v : J — > H such that 
the restriction v\i n extends to an r-times continuously differentiable function on the 
closed interval [i n _i,i„], for 1 < n < N. In other words, v is a piecewise C r function 
with respect to the time levels t n . 

If v € C 1 ( J, H), then its fractional derivative (jl.2l) admits the representation [31] 

™-l « /.min(tj,t) 

B a w(t) = w 1+Q (i)v° +Y / oj 1+a {t-t J )[vY / w 1+a (t - d« (2.3) 

for i G /„ and — 1 < a < 0. Thus, B a v(t) is left-continuous at t = t n _i but has a 
weak singularity (t — i„-i) a as t — > if [u]™ -1 ^ 0. However, if < a < 1 then 
B a v(t) is continuous for < t < T. For — 1 < a < 1, the piecewise-linear DG time 
stepping procedure determines U G W by setting = ito and requiring (301 131] 

(U»- \ Xl~ + jf [(I7'(t), X(t)) + A(S a [/(t), X(t))] <K 

= (CP- 1 , A^ 1 ) + jf </(i), X(t)> dt, (2.4) 

for 1 < n < iV and for every test function I 6 Pi. The nonlocal nature of the 
operator B a means that at each time step we must compute a sum involving all 
previous times levels, but this sum can be evaluated via a fast algorithm [2D] , 

2.3. Galerkin orthogonality and stability. For v G C 1 ( J, £>(A 1 / 2 )) and to G 
C(J, -©(A 1 / 2 )), we define the global bilinear form 

N-l n . 

g n (v, w) = (v° + ,w° + ) + ([<> o + E / w > + w )] ( 2 - 5 ) 

Summing the equations (|2.4[) gives 

G N {U,X) = (U°,Xl)+ / (f(t),X(t))dt foraUXGW, (2.6) 

Jo 

and conversely, (|2 .6[) implies that {J satisfies (|2.4[) for 1 < n < N. Since [w]™ = 0, 

G N {u,X) = (u ,Xl}+ (f(t),X(t))dt, (2.7) 
Jo 
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and thus, assuming U°_ = u , the error has the Galerkin orthogonality property 

G N (U-U,X) = forallXeW. (2.8) 

The DG method is unconditionally stable. Indeed, with the notation 

||[/|U = sup||C/(i)[| for any JC [0,T], 
tei 

the following estimate holds. 

Theorem 2.1. Given U° e H and f e Li((0,T);H), there exists a unique 
U e W satisfying (|2~4l) for n = 1, 2, N. Furthermore, U(t) £ D(A) for t > 0, 
and 



\U\\j n < 



(WL,K)+ / W),U(t))dt 



forl<n<N. 



Proof. Since (V',V) = ±(d/dt)\\V\\ 2 and A(B a V, V) dt > we find that 



N-l 



2G N {V,V) > \\V°\\ 2 + \\V N \\ 2 + \\[ V } n f for all yeW, (2.9) 



n=l 



(2.10) 



implying the stated estimate [30l Theorem 2.1], [3ll Theorem 1]. □ 
2.4. A discontinuous quasi-interpolant. The conditions 

U~v(t-) =v{t~) and J [v(t)-ITv(t)]dt = 
determine a unique projection operator II - : C(J, H) — > W. Explicitly, 
U-v{t) := v(t~) + Vn ! (t - t n ) for t £ I n , 

K n f Z 

where v n = k~ x L v(t) dt denotes the mean value of v over /„, and the interpolation 
error admits the integral representations [30, Equation (3.8)] 



n~«(t) - v(t) = I v '(s)ds~2\^ 

Ki j i 



(s - t n -i)v'(s) ds 



(t - s)v" (s) ds + 



tn t 



1.2 



n J I 



(s — t n _i) v"(s) ds, for t £ I n . 



Likewise, the conditions 

H + «(t+_i) = Ktf-x) and jf [v(t)-ll + v(t)}dt = 0, 
determine a unique projector n + : C(J, H) — > W, with 



(2.11) 
(2.12) 



k n /2 



{t - tn-l) for t e I n 
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and 

IL + v(t) - v(t) = - f v'{s)ds + 2 t ~^ 1 I (t n -s)v'(s)ds 



(s-t n _ 1 )v"(s)ds+ t —^ I (t n -s) 2 v"(s)ds. 



(2.13) 



Thus, short calculations lead to the error bound 
lin^-wHi, < (4-r)^- 1 / \\v {r) (t)\\dt < (4-r)kl\\v^\\ In far re {1,2}, (2.14) 



and the stability estimates 

lin^ik^Hk, iicn^^H^^Ajf \\v'(t)\\dt, Hp^m^ \\v'(t)\\dt 

(2-15) 

3. Dual problem. 

3.1. Properties of the adjoint operator. The adjoint operator appearing in 
the dual problem (| 1 . 5|) should satisfy, for appropriate u and v, the identity 

{v,B a w)dt= [ (B*v,w)dt, (3.1) 
Jo 

and the next lemma establishes an explicit representation of £>* . 
Lemma 3.1. The identity (|3.ip holds in the following cases. 
1. If —1 < a < and u, w e C^J, H), with 

d f T 

8*w(t) = — — / wi+ Q (s — t)w(s) ds fortEJ. 



dt,. 

2. J/0 < a < 1 and u, w e C(J, H), u»ift 



B* a w(t)= I u a (s ~ t)w(s) ds forO<t<T. 



Proof. In case 1, we see from the representation (|2.3[) that 

/ (B a v, w) dt = y~] / (B a v, w) dt = Si + S 2 + S3 
Jo n=1 Ji„ 

where, letting J?"' J = Jj uji +a (t — tj)w(t) dt and 

/> />min(ij,i) 

D nj = / / uj 1+a (t - s)(v'(s),w(t)) dsdt, 



we dehne 

iV N n-1 N n 

n— 1 n—2j—l n—lj — 1 
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By reversing the order of integration, integrating by parts and then interchanging 
variables, we find that for 1 < j < n — 1, 

D nj = (vl,B n >i) - (iP-\B n j- 1 )- jf (v(t),^J^ LJ 1+a (s-t)w(s)ds^dt, 

whereas D m = —(v"~ 1 ,B n ' n ~ 1 ) — j In (v(t), §- t J* n w 1+a (s -t)w(s) ds) dt. Thus, after 
interchanging the order of summation for the double integrals, 

N N n-1 N 

n=l n =2 j=l 3 = 1 j 

that is, 

jV n-1 AT n— 1 ~T 

n=2 3=1 n=2 j=0 ^° 

= -S 1 -S 2 - [ {v,B* a w)dt, 
Jo 

so (13. lp holds. In the case < a < 1, we simply reverse the order of integration. □ 
The adjoint operator admits a representation analogous to (12.31) . 
Lemma 3.2. If — 1 < a < 0, then B* a w(t) equals 

N-1 N tj 

U>l+a(t n ~ t)w N - ^ Wl+ a (tj ~ t)[u>] J ~ ^ / Wi +Q (s - t)w' (s) ds 

j=n j=n*' niax (*i-i>*) 

for w £ C ll (J, H) and £ € /„. Thus, B^w(t) is right- continuous att = t n but possesses 
a weak singularity (t n — t) a as t —> t~ . 

Proof. li t £ I n and j > n + 1, then, integrating by parts, 

Wl+ a (s-t)w(s) ds = L0 2 +a{t ] -t)w j _ -UJ2+a(tj-l~t)w j + 1 - j U)2+a (s - t)w'(s) ds 



and J. ra u)\-\- a (s— t)w(s) ds — uj2+ a (t n — t)u>™ — f^ u;2+ Q (s — t)w'(s) ds. Differentiating 
these expressions with respect to t, we see from part 1 of Lemma l3.1l that B^w{t) equals 

N N N . tj 

^u>i +a (tj —t)wP_ - ^2 vx+ a {tj-i -t)w^ 1 / u) 1+a (s-t)w'(s)ds, 

j=n j=n+l 3= „Jmax(t i _ 1 ,t) 

and the result follows after shifting the index in the second sum. □ 

3.2. Representation of the nodal error. Integration by parts in (|2.5I) . to- 
gether with the identity (I3.1[) . shows that for all v, w 6 C 1 (J, H), 

N-1 N 

N ...N 



^2(v!,[w] n ) + 1t [ [-(v,u')+A(v,B* a w)]dt. (3.2) 

71=1 11=1 "'•'•» 



Since — (v, z') + A(v,B^z) = (v,—z' + B^Az) = 0, the solution z of the dual prob- 
lem (|1.5p satisfies 

G N (v,z) = (v N ,z T ) for all v e C{J,D{A 1 ' 2 )). (3.3) 
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We therefore define the DG solution Z G W of (11.51) by 

G N (V, Z) = (V N , Zf) for all V eW, (3.4) 

with = zt, and deduce the Galerkin orthogonality property 

G N {V, Z-z)=0 for all V eW. (3.5) 

The following representation is the basis for our analysis of the nodal error. 

Theorem 3.3. If u and z are the solutions of the initial-value problem (|l.ip and 
of the dual problem (|1.5p . and if U and Z are the corresponding DG solutions, then 

(E7_ — u{tji), zt) — Gn(u — ILu, Z — z) for every zt G H. 

Proof. Taking V = U in (|3.4j) and v — u in (|3.3[) gives 

{U N - u(t N ), z T ) = (U N , z T ) - (u(t N ), z T ) = G N {U, Z) - G N {u, z) = G N {u, Z - z), 

where the last step used the Galerkin orthogonality property (|2.8|) of U, with X = Z. 
Now use the Galerkin orthogonality property (|3.5t of Z, with V = Il~u. □ 

3.3. Error in the DG solution of the dual problem. We will use the fol- 
lowing regularity estimates. 

Lemma 3.4. For —1 < a < 1 and < t < T, the solution z of the dual 
problem (|1.5j) satisfies 

IIA-V(t)n + (t — tjiiA-vwii < c(t - *r imi 

and 

(T - + INWII + (T- t)\\z'(t)\\ < C\\z T \\. 

Proof. Define the time reversal operator 1Zv(t) — v(T — t). Since IZdt — —d{R, 
and 1ZB* = B a lZ, we deduce from (11.51) that the function v = 1ZA~ x z satisfies 

v' + B a Av = for < t < T, with v(0) = A~ x z T - 

Known results for —1 < a < [19l Theorem 4.2] and < a < 1 [21j Theorem 2.1] give 
\\v'(t)\\ +t\\v"(t)\\ < Ct a \\Av{0)\\ = Ct a \\z T \\, implying the first estimate. Similarly 
[T^l Theorem 4.1], the function w = IZz satisfies 

t 1+a \\Aw(t)\\ + \\ w (t)\\+t\\w'(t)\\ < C\\w(0)\\ = C\\z T \\, 

implying the second estimate. □ 

To investigate the DG error for the dual problem, we make the splitting 

A- 1 (Z-z) = C + 6 where C = - z) and 6 = A^ 1 (Z - LT + z) G W. (3.6) 

Lemma 3.5. The function ( in (J3S]) satisfies < Ct% + k 1+a - \\z T \\. 
Proof. By (|2.14|) and Lemma EOI ||C||j is bounded by 

||(Z-n + )A" 1 z|| 7 < 3 max / 1 z'(t)\\ dt < C\\zt\\ max / (t N ~t) a dt. 

l<n<NJ In l<n<NJ In 
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If -1< a < 0, then (1 + a) J In (t N - t) a dt = (t N - t n -\) l+a - (t N - t n ) 1+a < k}+ a , 
whereas if < a < 1, then J. (tpf — t) a dt < k n t%. □ 

Lemma 3.6. The function 9 in (EU) satisfies ||0|| 2 7 < 8\J^ (Q(t), B* a AC,) dt\. 
Proof. By ([33|), G N {V, C + 9) = GiA^V, Z - z) = for all V € W, where we 
used the identity Gjv(u, j4 _1 k;) = Gjv(j4 _1 v, u>) and the fact that A~ l V e W. Thus, 

G JV (V r , 9) = -G(V, C) for all V e W. 

Since = for < n < N - 1, the formula (02]) shows 

AT AT „ 

n=l n=l 7 ™ 

and integration by parts gives J In (V,£')dt = (VT,C-> ~ fi n ( V ',O dt = (V-X-), 
where, in the last step, we used the second property in (|2.12l) and the fact that V is 
constant on I n . Thus, if we define g = — AB^C = —B„A( then 

G N (V, 6) = - f A{V, B* a C) dt= [ (V, g) dt for all V eW, 
Jo Jo 

which means that 9 € W is the DG solution of -& + B*A8 = g(t) for < t < T, 
with 6(T) = 0. The desired estimate follows by the stability of 9, which we can prove 
by applying Theorem O to 72.6 (t) = 6(T - t). □ 

Recall that t(t) = max(l, | logt|). 

Lemma 3.7. If — 1 < a < 1 then 

T 



(V,B* a AC)dt 



< Cr N + k L+a -£(t N /k N )\\V\\j\\z T \\ for all V e W. 



Proof. Suppose first that — 1 < a < 0. Since B a V — {B\ +a V)' and C™ 1 = 0, we see 
using p.l[) and integrating by parts that 



N N 



f (V,B* a A()dt = f2 [ ((Bi +a Vy,A()dt = f2 [ (A n , AC') dt, 
Jo n=1 Ji n n=1 Ji n 

where, for i£l„, 

A"(t)=B 1+a F(t„)- J B 1+a y(t) 

= / [w 1+a (* n - s) - - s)]V(s)ds + / Wi +a (tn - s)V(s)ds. 

Jo Jt 
The function wi+ Q is monotone decreasing whereas uJ2+ a is monotone increasing, so 

||A n (t)|| < \\ v h(J \ui+a(t n -a)-ui+ a (t-8)\da + J UJ 1+a (t n ~ s) ds^ 

= \\V\\j[uJ2 +a (t)~U2 +a (t n )+2u2 +a (tn-t)] < 2 1| V|| jU) 2 + a (t„ - t) 

and the Cauchy-Schwarz inequality shows that |/q T (V, £>* -AC) is bounded by 
2\\V\\j(^2^+ a (kn) \\AC\\dt + J^ w 2+a (t N -t)\\A('\\dt 
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The integral representation of the interpolation error (|2.13p and Lemma 13.41 imply 

JV-l „ N-l 



X; [ \\ac\\ dt<cJ2 k n +a I \\A 



(It 



T-k? 



<Cfc i+tt ||z T || / (T -t)- 1 dt = Ck 1+a log p- 
Jo k N 

and f lN io 2+a {t N -t)\\AC\\dt < C\\z T \\ J In (t N - t) a dt < C\\z T \\k]+ a . The desired 
estimate follows at once. 

Now let < a < 1. By part 2 of Lemma \'3.1l 



T 



(V,B* a A()dt 



< / \\V(t)\\ / u a (s-t)\\A{(s)\\dsdt 



<\\V\\j \\A({s)\\ / uj a (s -t)dtds = \\V\\j I \\AC(s)\\u 1+a (s)ds 
Jo Jo 

<CT a \\V\\ J [ umidt. 



The estimates (I2.14[) and (|2.15j) imply that 

~t n n-1 r 

/ \\Aat)\\dt<^k n \\AC\\ In <4k N \\z\\ lN +3j2k n / WAM dt, 

Ja n=l n=l ^" 

and we know from Lemma 13.41 that < C||zt|| and 

N—l . ptN-l 

Vfc„/ \\z'{t)\\dt<Ck n \\z T \\ {t N -t)- 1 dt = Ck N \\z T \\\og{t N /k N ). □ 

n=l ^0 

Hence, we arrive at the following error estimate for the dual problem. 
Theorem 3.8. Let z denote the solution of the dual problem (jl.5l) . and let Z 
denote the DG solution defined by (|3 .4[) . Then, for — 1 < a < 1, 



A-\Z - z)\\j < Ct a N + k 1+a -e(t N /k N )\\z T \ 



Proof. The splitting J3I5J implies that \\A- 1 {Z - z)\\j < + \\Q\\j, and we 
estimate these two terms using Lemmas 13.51 13.61 and 13.71 □ 

4. Nodal superconvergence. With the help of Theorems 13.31 and 13.81 we are 

now able to estimate the error in the approximation £7™ rs u(t n ). Define 

n 

e{u) =D 1 +E 1 + max kjDj + V k^-Ej, (4.1) 

where 

Di = / ||Au'(i)||d< and E x = [ t 1+a - \\A 2 u' (t)\\ dt, (4.2) 

Jlx Jlx 
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with 



D n = 



\\Au"{t)\\dt and E n=j i ll^ 2 ""WII 



dt for 2 < n < N. 



(4.3) 



Theorem 4.1. Let u be the solution of the initial value problem and let U 

be the DG solution satisfying (|2.4|l . Then, for 1 < n < N , 



\\U n -u(t n )\\ < Ct 2 n a +k i+a -£(t n /k n )e(u) 



(4.4) 



Proof. Put rj = u — II u and define 

8*1 = J (rj, {z - Z)') dt and 8% = y (B tt A?7, z - Z) dt. 
Since 77™ = for 1 < n < N, we see from Theorem 13.31 Lemma T3. II and (|3.2j) that 

N 

(U N - u(i n ), z T ) = G N ( V , Z-z) = Y,(&x+ 82)- 

71=1 

Since Z' is constant on I n , the second property of II - in (|2.10p gives 

8[ l = f (r),z')dt= f (n(t),z'(t)-z'(t n - 1 ))dt= f f (Ar ] (t),A- 1 z"(s))dsdt, 

Jin Jin Jin Jtn-1 

and therefore, using Lemma |3.4[ 

N-l N-l . N-l 

£ I*? I < II^IU £ *» / \\A- 1 z"(t)\\dt<c\\Ar ) \\j N \\z T \\'£k n / (t N -t) a - l dt, 

71=1 71=1 71=1 

whereas |<*f | = |/ 7 (A77, A~ 1 z') dt\ < C\\Ar]\\ lN \\z t \\ fj {t N ~t) a dt. Here, 



71=1 



J2kn / (t N - t)*- 1 dt < k (tN-t^dt^-it^-k^^Ct^+k 1 ^-, 



and likewise J /iv (t A r-i) Q rfi = /c] v +Q /(l + a) <C^+fc 1+a -. By (|2~14D . ||A?7|| Zl < 3£> a 
and mry]]z n < 2k n D n for 2 < n < N, so 

£|*TI <c4 Q+ /c 1+Q - ||*r||(l>i+ max k n D n \ for -1< a < 1. (4.5) 

n=l ^ 2 _ n - y 



Turning to <5f , if -1 < a < 0, then [3U Lemma 2] 



1V ft N / jv 

£<5£ = / (B a A 2 ri, A~ r (Z - z)) dt < C\\A~\Z - + 

n=l ^ n=2 

but if < a < 1 then j^ 1 1 is bounded by 



\\A-\Z-z)\\ Iri \\B a A 2 r,(t)\\dt< \\A-\Z-z)\\ In / / w a (t-s)||AVs)||rfs, 
J/,, J/„ Jo 
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so, after summing over n and reversing the order of integration, 



Y d m<Ct%\\A-\Z-z)\\ J ( 

n=l J ° 

The integral representation (|2.11j) implies that 



A 2 ri{t)\\dt. 



\\A 2 v{t)\\< I ( f 1 \\A 2 u'( S )\\ds + 2%^ I s\\A 2 u'{s)\\ds)dt 



h -Hi 



I \\A z u'(s)\\ds + 2^ vr - / s\\A z u\ 
t fc i Jii 

J \\A 2 u'{s)\\^ dt + ^J^ 2(ti - t) dtj ds = 2E 



and by (EH, fi» ||A 2 »,(i)ll* < £n= 2 M^lk < En=2 2k l E J- Applying Theo- 
rem 13.81 

W I ^ Ctf + k 1+a -£(t N /k N )\\z T \\ (e x + k 2 + a -E n ) for -1< a < 1. (4.6) 

n=l ^ n=2 ' 

Since zt S H is arbitrary, the desired estimate follows from (|4.5p and (|4.6p . □ 

To estimate the convergence rate at the nodes, we introduce some assumptions 
about the behaviour of the time steps, namely that, for some fixed 7 > 1, 

K < C 7 fcmin(l,4~ 1/7 ) and t n < C 7 i„-i for 2 <n < N, (4.7) 

with 

c 7 fc 7 < fei < C 7 fc 7 . (4.8) 
For example, these assumptions are satisfied if we put 

t n = (n/iV) 7 T for < n < N. (4.9) 



Lemma 4.2. Assume that u satisfies (|1.3|) and (|1.4j) . and that the time mesh 
satisfies (|4T71) and ([ITS]) . TTien, mi/i 7* = (2 + <7 and /or 1 < n < AT, 

'fc 7CT , 1<7<7*, 
<(</>< C T M x <( k 2+a -l{t n /k x ), 7 = 7*, 
fc 2+Q -, 7>7*- 

Proof. The stated assumptions imply that D x + Ex < CM J* 1 t"- 1 dt < CMkf < 
CMW, and, for 2 < j < n, 

f [ k la \ < 1 < 21 a 

kjDj < CMkj / t a ~ 2 dt < CMk 2 t°- 2 <CMxl „ ~ 1 1 ' 

r fc 7CT , 1 < 7 < 7 *, 

<CMxlk 2+a ~, 7*<7<2/ct, 

l < jr-(2+a_)/ 7jfc a +a _ j 7 > 7 ^ 
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Similarly, 



3=2 3=2 



2+ct_ 



x(T— 3— a- 



dt < CMk 2+a - 



t a-l-(2 +a -) h dt 



< CMk 



2+a_ 



' fc7 <r-(2+a_) ; 1< 7<7 * J 

< < log(f„/fei), 7 = 7*, 



We can now state our main result on nodal superconvergence. 

Theorem 4.3. Assume that the solution u of the initial value problem (|1.1[) 
satisfies (|1.3j) and fj 1 . 4 [) . and t/iai i/ie iirae mes/i satisfies (|4.7|) and (|4.8|) wraf/i 7 > 
7* = (2 + a_)/a\ Then, for the DG method (|2.4[) . we /lave i/ie error bound 



\\U1 - u{t n )\\ < CMt 2 n a +k 3+2a -e(t n /k n ) forl<n<N. 



Proof. The error bound follows at once from Theorem 14.11 and Lemma 14.21 □ 

5. Postprocessing. We can postprocess the DG solution U to obtain a globally 
superconvergent solution {/" using simple Lagrange interpolation, as follows. Given a 
piecewise continuous function v : J — > H, define Cv : J — > H by linear interpolation 
on the first two subintervals, 

Cv(t) =k- 1 [{t n ~t)v n r 1 + (t-t r ^ 1 )v n :} for t e I n and n £ {1,2}, (5.1) 

and backward quadratic interpolation on the remaining subintervals, 

£ V U\ _ ft - t n -i)(t - t n ) ^„_ 2 _ {t - t n -2)(t - t n ) ^ n _ x ^ (t - t n - 2 )(t - t n -i) ^ n 

kn—x{kn—\ ~t~ ^n) k n —\k n {k n — 1 ~\~ k n )k n 

(5.2) 

for t 6 7„ and n > 3. Thus, (£i;)(i n ) = u" for < n < N, and we define the 
postprocessed solution by 

[/» = CU. (5.3) 

The interpolant of the exact solution satisfies the following error bound. 
Lemma 5.1. If there exist positive constants M and eft such that 

\\u'(t)\\+t 2 \\u"'(t)\\<Mt a *- 1 forO<t<T, (5.4) 

and if the time mesh satisfies (|4.7[) and (j4.8[) with 7 > 3/o"", then 



W a \ 1<7<3/o- b 



||u-£u||j < CM x 



Proof. If n e {1, 2} and t € I n , then 



t 4 ,t„ 



(u — Cu){t) = u'(s)ds — — / w'(s) ds 
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and thus ||u - Cu\\i n < 2 f*"^ \\u'(s)\\ds < CMt°* < CM(ki + k 2 y" < CMP 1 , 
If n > 3 and t £ I n , then we can write the interpolation error in terms of a divided 
difference, (u - Cu)(t) = u[t n -2,t„-i,t,t n ](t - t n - 2 )(t - t n -i)(t - t n ), so 

||u - Cu\\ In < \k 2 n (k n ^ + fcn)^||u w || [tB _ ajtn] < CMk 2 n {k n ^ + k n K*~ 3 , 

where, in the final step, we used (|4.7[) . If 1 < 7 < then, again using (14. 7} , 

k 2 n (k n ^ + kk*- 3 < cikt^r^k^t**- 3 = Ck^(k n /t n ) 3 -^ < Ck^\ 

but for 7 > 3/cr, 

k 2 n {k n . x + k n )t<^~ 3 < Cikt^hft^- 3 < Ck 3 t^- 3 h < CT° f - 3 ^k 3 . □ 

Now consider the stability of the interpolation operator C. We see from (|5.1j) that 

\\Cv\\ h <max(\v°_\,\vh\) and ||£«||j 2 < max(|ui|, \v^\). 

A similar estimate holds for the subsequent subintervals provided the mesh satisfies 
the local quasi-uniformity condition 

k n < Afc„_i for 3 < n < N. (5.5) 



For example, our standard mesh (|4.9p satisfies this condition with A = 2 7 — 1. 
Lemma 5.2. If (|5.5|) holds, then 



\\£v\\i n < (2+ |A) max \vt\ for3<n<N. 

71 n — 2<j<n 



Proof. The estimate follows from (|5.2j) because, for t G I n and n > 2, 

\(t - t n -l)(t - t n )\ < < jfcn < 1 . 

k„-i{k n -i + k„) ~ fc„_i(fc„_i + fc„) ~ k„-i ~ 4 

|(t - t n - 2 ){t - t n )\ < (fc ra _i + fc»)fcw _ ]_ , fc » < 1 + A 
\(t - t n - 2 )(t - t K _i)| < (fcn-i + fc„)fc„ _ _^ 

(k n — 1 + kn)k n (kn—i + k n )k n 

□ 

Hence, the interpolant ?7* is superconvergent, uniformly in i. 
Theorem 5.3. Suppose that the time mesh satisfies (|4.7|) . ()4. 8|) and (|5.5|) . and 
that u satisfies (j5.4[) . J/ 7 > 3/a" , i/iera i/ie postprocessed solution (15 ,3|) satisfies 

||f/» - u[[j < max ||tT - u(t„)|| + C T AM x ' 1 ; ^j 3 ^' 

0<n<N I k°, 7 > 3/C7 B . 



Proof. Write [/» - u = (C u — u) + £(£/ — u) and apply Lemmas 15.11 and 15.21 □ 
6. Spatial discretization. 
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6.1. The fully discrete DG method. We denote the norm of u in H r (£l) 
by \\u\\ r , and assume now that A = — V 2 in a bounded, convex or C 2 domain f2 in M. d , 
subject to homogeneous Dirichlet boundary conditions. Thus, if u G Hq(£1) and 
Au G L 2 (Q), then u G ff 2 (ft) and ||w|| 2 < C[[Ait[[ ia(n ). Let C D{A 1 / 2 ) = H^(ft) 
denote the space of continuous, piecewise-linear functions with respect to a quasi- 
uniform partition of f2 into triangular or quadrilateral (or tetrahedral etc.) finite 
elements, with maximum diameter h. Recall that the -^-projector Ph : L 2 {^1) — > Sh 
and the Ritz projector Rh : Hq(&) —> Sh are defined by 

{P h v, W) = (v, W) and A(R h v, W) = A(v, W) for all W G S h , (6.1) 

and that the latter has the quasi-optimal approximation property 

\\v-R h v\\+h\\V(v-R h v)\\<Ch 2 \\v\\ 2 for v G Hq(Q) R H 2 (fl). (6.2) 

Let W(Sh) denote the space of piecewise linear functions U : J — > Sh (so U is 
continuous in space, but may be discontinuous in time). 

We define the fully discrete DG solution Uh G W(Sh) by requiring (12.41) to hold 
for every X G W(S h ). Equivalently, cf. (|2Hj) . 

G N (U h ,X) = {U%_,X°) + (f(t),X(t))dt for all X G W(S h ), (6.3) 

Jo 

where, for simplicity, we choose J7°_ = (Uh) . = PhUo- In view of (|2.7p . the Galcrkin 
orthogonality property (|2.8p now takes the form 

G N (U h -u,X) = {P h u -u ,X + ) =0 for aU X G W(S h ). (6.4) 

Similarly, the fully discrete DG solution Zh G W(Sh) for the dual problem (| 1 . 5[) is 
defined by 

G N (V,Z h ) = {V N ,Z^ + ) for all V G W(Sh), with Z% + = P h z T , (6.5) 

and, since z satisfies (|3.3|) . 

Giv(V, Z h - z) = (Vj^, PhZT - z T ) - for all V G W(S h ). (6.6) 

Theorem 13.31 generalizes as follows. 

Theorem 6.1. If u and z are the solutions of the initial value problem (jTTTJ) 
and the dual problem (|1.5[) . and if Uh and Zh are the corresponding fully discrete DG 
solutions satisfying (|6.3p and (|6.5p . i/ien 

(f7/f_ — u(tw), 2t) — Gn(u — T\T RhU, Zh — z) for every zt G H. 

Proof. By taking V — Uh in (|6.5p we find that 

(C/f_,z T ) = (U?_,P h z T ) = (U^Zb) = G N (U h ,Z h ), 

and taking w = u in (|3.3p we have {u{tjq), zt) = Gn{u, z), so 

(U^- - u(t N ),z T ) = G N (U h , Z h ) - G N (u,z) 

= G N (U h - u, Z h ) + G N (u, Z h - z) = G N (u, Z h - z), 

where the final step used (|6.4p with X = Zh- Since 

G N (u, Zh—z) = G N (u - Tl~R h u, Z h — z) + G N (U^R h u, Z h - z), 

the result follows after putting V = H~RhU in (|6.6p . □ 
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6.2. Error in the fully discrete DG solution of the dual problem. We 

modify the splitting (|3.6I) . by writing A~ 1 (Zf l — z) = ( + ^ + ^ where 

C = A~ 1 (U + z — z), = A- 1 U+(P h z-z), <P = A- 1 (Z h -n+P h z)eW(S h ). 



Theorem 13.81 generalizes as follows. 

Theorem 6.2. Let z denote the solution of the dual problem (|1.5I) . and let Zh 
denote the fully discrete DG solution defined by (|6.5p . Then, for —1 < a < 1, 

\\A-\Z h - z)\\j < c(t a N + k 1+a - + h 2 )e(t N /k N )\\z T \\. 



Proof. We already estimated ||£|| in Lemma [3.51 To estimate observe that 
since A -1 commutes with LT + and since ARh = PhA (implying A~ l Ph = RhA^ 1 ), 

* = U + {R h -l)A- 1 z. (6.7) 

Using (|2.15p . the error bound (|6.2I) for the Ritz projection, iJ 2 -regularity for A and 
Lemma EOl we find that 



||*||/„ < 3||(i2fc -VjA^zW^ < Ch 2 \\A- l z{t)\\ 2 < Ch 2 \\z(t)\\ < Ch 2 \\z T \\. (6.8) 
To estimate observe that since A~ Y V = A~ x P h V = RhA~ l V, 
G N (V, C + * + $) = GniA-ty, Z h -z) = G N {R h A- l V, Z h -z) = 0, (6.9) 
where we used (|6.6[) with V replaced by R/ l A~ 1 V. From the proof of Lemma [3T6l 

G N (V,0 = [ T (V 7 B* a A0dt, (6.10) 
Jo 

and by ([3~2")l . 

N-l N . 

G N (V,y) = (V N ,* N ) + Y,(V-,m n ) + Y, [-(V,*')+A(y,B* a *)]dt. (6.11) 
Since A(V, B*j!) = A(V, 3*^X1+ {P h - X)z) = {V, (P h - l)B*Jl+z) = 0, 

\G N (y,n< ii^iui 11^11 + E ii[*ni +E / n*'iH- 

^ n=l n=l"' / ™ / 

By dUHl), ||*_ || < < Ch 2 \\z T \\. Using (|67f| . (I2TT51) . (IQl and LemmaGEl we 

have ||[tf] n || < / Jb [[(iZ fc -2)A-V(t)[|dt < C7i 2 ||z r || / Jn (T - t)" 1 dt so 

^-1 /-tjv-i 

J2 \M n \\ <Ch 2 \\z T \\ / (T-t)- 1 dt = C/ l 2 ||^ T ||log(i JV Aiv). 

n=l 

Using ((2~T5|) . ([67T]1 and Lemma El we find that 

/ W(t)\\dt<^f (t n -t)\\(R h -l)A- l z'(t)\\dt<^\\z T \\ [ %r^dt, 

Jl n K n J I n K n Jl n 1 — I 
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so 

f f f^N—i r \ 

Y J J i W\\dt<Ch 2 \\z T \\^J^ {T-t)- x dt + kx X j i dtj, 

and we conclude that |Gjy(V,*)| < Ch 2 £(t N /k N )\\V\\j\\z T \\. Therefore, by (j6~9)) . 
(j6~TUj) and Lemma 

\G N (y,*)\ = \G N (V,Q+G N (V,*)\ < C(t a N + k 1+a -+h 2 )£(t N /k N )\\V\\j\\z T \\. (6.12) 
Fix n with 1 < n < N, and define V G W(S'/ l ) by 




0, if t € ij for 1 < j < n- 1, 
if t G for n<j<N. 



In view of (HU, G^V,*) > i||$^|| 2 + UK'T + k £f=n WP. 80 the esti " 
mate {6T2) gives H^ 1 ]] 2 + ||[$]"|| 2 < C{t% + k l + a ~ + h 2 )£(t N /k N )\\$\\ [tn _ uT) \\z T \\ 
for 1 < n < N — 1, whereas 

piV-Y + ||# jv ||2 < c-( t ^ fc i+«- + /l 2 )i( tjv /A ;w )||$|| (tN _ 1)T) ||z T ||. 
Furthermore, 1 1 <I? 1 1 = max (||$+ _1 ||, II*™ ||) because $ is piecewise linear in t, and 

mw < wnw + win implying that n$n 2 n < i^r 1 !! 2 + imii 2 + pni 2 - % 

letting n* = argmax 1<n<w || $||/„, we see that 

ii* ii 2 / = 11*111. < c(t a N +k i+a - + h 2 )i{t N /k N )\\n j\\zti 

giving the desired bound for ||$||j. □ 

6.3. Fully-discrete nodal error. As claimed in the Introduction, we have the 
following error bound for Uh- 

Theorem 6.3. Assume that the solution u of the initial value problem (|1.1[) 
satisfies (jl.3[) and (|1.4[) . and that the time mesh satisfies assumptions (14.71) and (I4.8P 

with 7 > 7* = (2 + a_)/er. Then, the fully discrete DG solution Uh G W(Sh) satisfies 

\\U%_ ~ u(t n )\\ < C T M(k 3+2a -l(t n /k„) + h 2 ) forO<n<N. 

Proof. In view of Lemma 14.21 it suffices to show (cf. Theorem 14. ip that 

||L7£_ - u(t n )\\ < C T {k 1+a - + h 2 ) £(t n /k n )e{u) + C T Mh 2 . 

Put £ = u — RhU and rj = u — Tl~u so that u — II RhU = rj + LT~£ and thus, 
by Theorem O (Ug_ - u{t N ),z T ) = G N { Vl Z h - z) + G N (J[-^Z H - z). Using 
Theorem 16.21 in place of Theorem 13.81 we can show as in the proof of Theorem 14.11 
that \G N ( V ,Z h -z)\ < C T \\z T \\ (k l+a - + h 2 )£(t n /k n ) e{u). By & G N (IL-£,Z h ) 
equals 

N-l n . 

((n-o°+, z° h+ ) + J2 (p-fl", z% + ) + J2 K(n^)', z h ) + A(B a ir& z h )\ dt, 

n=l n=l"' / " 



18 



KASSEM MUSTAPHA AND WILLIAM MCLEAN 



and, since B a commutes with the Ritz projector , the definition (I6.1[) of implies 
that A(23 Q II-£, Z h ) = A(B a n-u, Z h ) - A(R h B a nru, Z h ) = 0. Integrating by parts, 
applying the interpolation and orthogonality properties (12.101) of II - , and noting that 
£™ = i{t n ) = (n _ £)™ and that Z' h is constant on /„, 

jf ((u-a,z h )dt = ((u-o^,z^)-((u-o n + - 1 ^i:+ 1 )- 1 (n-£,z' h )dt 

= (er 1 - (n-or 1 .^ 1 ) + / <c, ^> * 

so Gjv (n — z h ) = <£(o), + E^i^K]". ^+) + Eli f In ^) using m 

with ?; = n _ £, and noting that [£]" = 0, we obtain 

iV 

G N (n-e,^-^) = (C(o),^+)-(e(T),z T ) + ^ / (^,Z h )tft. 

n=l 

Stability of the fully discrete dual problem, < C||zr||, follows from ()2.9j) . so 

\G N (nr£, z h - z )\< c\\z T \\ (\\t(p)\\ + U(t N )\\ + J <ft) 

< Ch 2 \\z T \\ (\\Auo\\ + \\Au(T)\\ + £ \\Au'(t)\\ dt) , 

where we used the error bound (|6.2I) for the Ritz projector. The result follows using 
the regularity assumption (|1.3I) . □ 

6.4. Postprocessing the fully discrete DG solution. Theorem 15.31 remains 
valid if — CU and Z7" are replaced by lj\ = CUh and U£_, respectively. 

7. Numerical results. We present a series of numerical tests using a model 
problem in one space dimension, of the form (jl.ip with 

Au = -u xx , ft = (0,1), [0,T] = [0,1], u (x)=x(l-x), / = 0, 

and homogeneous Dirichlet (absorbing) boundary conditions. These tests reveal faster 
than expected convergence when a < 0, and that our regularity assumptions are more 
restrictive than is needed in practice. We apply the fully discrete DG method defined 
in Section HTTl employing a time mesh of the form (|4.9j) . for various choices of the mesh 
grading parameter 7 > 1, and a uniform spatial mesh consisting of M subintervals, 
each of length h = l/M. We always choose M — [TV 3 / 2 ] so that h 2 w k 3 and hence 
the error from the time discretization dominates the spatial error. 

7.1. The exact solution. Separation of variables yields a series representation 

00 

u(x,t) = 8^2lj~ 3 sm(uj n x)E 1+a (-L}lt 1+a ) with u> n = (2n + 1)tt, (7.1) 

where the Mittag-Leffler function is given by E u (t) = J2^Lo^ P 7^(1 + vp). We can 
verify directly that u satisfies the regularity conditions 

t 1+a \\Au'(t)\\+t 2+a \\Au"{t)\\ < Mt a - X for < t < T, (7.2) 
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Table 7.1 

The left nodal error maxi<n<iv \Wh— ~ u (*")ll an d the rate of convergence when a = —0.3, 
for different mesh gradings 7. 



N 


7-1 


7 = 2 


7 = 3 


7 = 3.25 


20 
40 
80 
160 


2.01e-03 
8.61e-04 1.220 
3.90e-04 1.143 
2.21c-04 0.821 


l.OSe-04 

3.15e-05 1.780 
9.33e-06 1.758 
2.77c-06 1.753 


6.39e-05 
1.09e-05 2.546 
1.81e-06 2.596 
2.92e-07 2.632 


6.39e-05 
1.10e-05 2.535 
1.82e-06 2.595 
2.94c-07 2.632 



Table 7.2 

The right nodal error maxo< n < at— 1 \\U£ + — u(t n )\\ and the rate of convergence when ct = —0.3, 
for different mesh gradings 7. 



N 


7 = 1 


7 = 2 


7 = 3 


7 = 3.25 


20 
40 
80 
160 


4.74e-02 
3.05e-02 0.636 
1.89e-02 0.689 
1.16e-02 0.710 


6.03e-03 
2.26e-03 1.416 
8.51e-04 1.410 
3.21e-04 1.406 


1.63e-03 

4.18e-04 1.966 
1.06c-04 1.982 
2.66c-05 1.990 


1.52e-03 

3.91e-04 1.964 
9.89e-05 1.982 
2.49e-05 1.989 



with 



||u(t)|| 2 +t|K(*)||2 <M for0<t<T. 
In fact, by differentiating ()7.1|) . 



(7.3) 



dlu xx {x,t) = -S^w^sinM- E 1+a (-uj 2 n t 1+a ) forj G {1,2}, 



n=0 



so by Parseval's identity, 



||^(i)|| 2 = |l^ x (t)f = 32^ 



x.r\ —E l+a {-^ n t 1+a ) 



The Mittag-Leffler function satisfies [THJ Theorem 4.2] 



< Ct- (1+a ^- j Lo- 2 >* for j e {1, 2,3,.. .} and |^| < 1, (7.4) 



and taking /i = — e yields 



00 c I f«(i+ti)+«v 

(t^ a \\d{Au(t)\\) 2 < Ct 2 < 1+a ^ +2a Y,^r 2 < -J- lor I . , i. 



n=0 



4c 



Thus, the regularity condition (|7.2I) holds for a = (1 + e)(l + a) < 4(1 + a) and 



M = C{\ — e) x l 2 . In particular, putting e = gives the bound for t||it'(t)||2 in (|7.3p . 
and since \E 1+a (~ujlt 1+a )\ < C for all * > we also have |K*)||| < C^o^ 2 < 
00. However, u fails to satisfy the second regularity assumption (jl.4p used in our 
theoretical analysis. 
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Table 7.3 

The left nodal error maxi<n<iv ll^fe— — an< ^ the rate of convergence when a = +0.3, 

for different mesh gradings 7. 



N 


7=1 


7 = 1.5 


7 = 1.75 


7 = 2 


20 
40 
80 
160 


2.10e-04 
6.77e-05 1.632 
2.19e-05 1.636 
7.116-06 1.625 


2.08e-05 

3.61e-06 2.527 
6.43e-07 2.486 
1.17e-07 2.461 


1.21&-05 

1.616-06 2.904 
2.13e-07 2.917 
2.80e-08 2.930 


1.23e-05 
1.57e-06 2.966 
1.99e-07 2.983 
2.53e-08 2.972 



Table 7.4 

The right nodal error maxo< n < at— 1 | Wh+ — u(t n )\ \ and the rate of convergence when ct = +0.3, 
for different mesh gradings 7. 



N 


7=1 


7 = 1.5 


7 = 1.75 


20 
40 
80 
160 


3.265e-03 
1.536e-03 1.088 
6.726e-04 1.191 
2.851e-04 1.238 


8.548e-04 
2.165e-04 1.982 
5.432e-05 1.995 
1.361e-05 1.997 


9.207e-04 
2.338e-04 1.977 
5.873c-05 1.993 
1.472e-05 1.996 



7.2. Nodal errors. The numerical results described below suggest that 

max ||^ -u(t n )\\<Ch*+Cx{ k Z' 1 ^^( 3 + "-)/ CT ' (7.5) 
l<n<iV H h K nJ>> \k 3+a -, 7>(3 + a_)/cr. V ' 

Thus, the time discretization error appears to be 0{k 3+a -) for 7 > (3 + a_)/<r, 
compared to our theoretical bound of 0(k 3+2a - ) for 7 > (2 + a-) /a, where the latter 
assumes the stronger regularity conditions (|1.3p and (|1.4[) . 

For a = —0.3, we observe in Table FTTI convergence of order / £ 1 - 25 7(a+i) f or \ < 
7 < (3 + a)/[1.25(a + 1)] ~ 3.086. In particular, the highest observed convergence 
rate is O(k 3+Ot - ), and not 0(k 3+2a - ) as expected from Theorem l6.3l Table I7T21 shows 
that the right-hand limit U£ + = U~h(t+) — Iim t _ > . t + Uh(t) is not a superconvergent 
approximation to u(t n ); the error is 0(k 2 ) at best. 

For a = +0.3, Table 17.31 shows convergence of order fc 125 7("+i) f or 1 < 7 < 
3/[1.25(a + 1)] w 1.85, so in the best case the error is 0(k 3 ), consistent with Theo- 
rem [Ol In Table [7T41 we see that U£ again fails to be superconvergent. 

Given a, it is natural to ask which value of 7 leads to the smallest error. Figure [7TT1 
shows the maximum nodal error (on a logarithmic scale) as a function of 7 G [1,8] 
for 4 choices of a, when M = 512 and N = 64 (so h 2 = k 3 ). The error is minimised 
when 7 w (3 + a_)/cr; for instance, in the case a = 0.2 the best choice is 7 ~ 
3/[|(1.2)] = 2. In Figure [7T2l we instead show the maximum nodal error as a function 
of a £ [—0.9, 0.9] for 4 choices of 7. The benefit from using non-uniform time steps is 
clear, except when a is close to —1 or 1. 

7.3. Global error after post-processing. We introduce a finer mesh 

g N ' m = { tj-! + ikj/m : j = 1, 2, . . . , N and £ = 0, 1, . . . , m}, (7.6) 

and define the discrete maximum norm ||u||j im = max t£ gN, m \\v(t)\\, so that, for suffi- 
ciently large values of m, ||£7 — u|| j jm approximates the global error ||J7 — u|| j. Now, in 
addition to the regularity assumptions (|7.2[) and (|7.3|) . we require that u satisfies 
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Fig. 7.1. The left nodal error maxQ<„<jv | \U£_ — 
-0.8, -0.4, 0.2 and 0.6, when M = 512 and N = 64 (so h 2 



u(t n ) 
= k 3 ). 



as a function of 7, for a 





Fig. 7.2. The left nodal error maxg< n <jv — u (tn)ll as a function of a, for 7 = 1, 2, 3, 

when M = 512 and N = 64 (so fe 3 = /i^ J. 




In fact, we see from (|7.1[) and (|7.4I) that, with fx = —1, 

00 

(^- 1 ||^)||) 2 < cr(t( 1 +-)- 1 ) a £ w ^ a < C^^- 1 ) 2 , 



so (pTl)) holds for cr* = 1 + a. Using Theorem [O] (cf. Subsection and (|73|) 
with cr" = (1 + a) < a « |(1 + a), we expect 

We observe this convergence behaviour in Tables 17.51 and 17.61 
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Table 7.5 

The uniform DG error after postprocessing, ||t/f — u\\ j t i%, and its rate of convergence, when a = 
—0.3, for different mesh gradings 7. 



N 


7=1 


7 = 2 


7 = 3 


7 = 3.9 


20 
40 
80 
160 


3.79e-02 
2.37e-02 0.675 
1.44e-02 0.716 
8.74e-03 0.724 


4.52e-03 
1.68e-03 1.425 
6.31e-04 1.416 
2.38e-04 1.410 


1.46e-03 

3.27e-04 2.154 
7.49e-05 2.127 
1.73e-05 2.113 


8.13e-04 
1.20e-04 2.763 
1.79e-05 2.743 
2.69e-06 2.735 


Table 7.6 

The uniform DG error after postprocessing, — m|| j,i2> and its rate of convergence, when a 
0.3, for different mesh gradings 7. 


N 


7 = 1 


7 = 1.5 


7 = 2 


7 = 2.35 


20 
40 
80 
160 


2.51c-03 
1.16e-03 1.120 
5.02e-04 1.205 
2.12e-04 1.245 


4.38e-04 
1.22e-04 1.845 
3.34e-05 1.867 
8.88e-06 1.911 


1.56e-04 
2.72e-05 2.515 
4.59e-06 2.568 
7.63e-07 2.588 


1.89e-04 
2.29e-05 3.046 
2.80e-06 3.029 
3.44e-07 3.024 



8. Concluding remarks. We have analysed a piecewise-linear DG method for 
the time discretization of (jl.ip — a fractional diffusion (— 1 < a < 0) or wave 
(0 < a < 1) equation — and proved superconvergence at the nodes, generalizing 
a known result for the classical heat equation. Numerical experiments indicate that 
our theoretical error bounds are sharp if a > 0, but not if a < 0. For generic regular 
data wo and /, derivatives of the exact solution are singular as t — > 0, but nevertheless 
by employing non-uniform time steps we achieve a high convergence rate of 0(k 3+a ~ ). 
After postprocessing the solution, the same high accuracy is achieved for all £, not 
just at the nodes. We have also proved that the additional error arising from a spatial 
discretization by continuous piecewise-linear finite elements is essentially 0(h 2 ). In 
future work, we aim to treat the case when the initial data uq is not smooth. 
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